Numerical integration using quad()

March 19, 2021

• First import some required modules

import numpy as np
from scipy.integrate import quad

An integral that comes up when trying to find the chemical potential of a system of identical Bosons is:

$\int_0^\infty \dfrac{x^2}{e^{x^2-\eta}-1}\, dx$

where $\eta=\mu/k_\mathrm{B}T$ and $\mu$ is the chemical potential.

This short tutorial deminstrate how to use quad() to numberically evaluate this integral for the $\eta = 0$ case.

• First, define a function that represents the integrand of the integral.

def integrand(x):
    return x**2/(np.exp(x**2) - 1)

• Next, we call the function inside quad() from the SciPy module to perfrom the numerical integration.

y, err = quad(integrand, 0, np.inf)

<ipython-input-27-472f14518a59>:2: RuntimeWarning: overflow encountered in exp
  return x**2/(np.exp(x**2) - 1)

• quad() returns two values. The first (y) is the value of the integral and the second (err) is an estimate of the error in the result.

y, err

(1.1575786866969464, 8.140518531748718e-09)

• Notice also that we got and overflow warning. This occurs because when $x$ is very large $e^{x^2}$ becomes exceedingly large. If we want, we can avoid this warning by simply truncating the integration using some maximum value of $x$, say $\infty\to 10$.

y, err = quad(integrand, 0, 10)
y, err

(1.1575786866970472, 7.989498536842354e-12)

For this choise of the upper integration limit, there was no warning and our result for y changed by only a completely negligible amount.